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A NIFMERICAL SIMULATION OF TRANSITION IN 
PLANE CHANNEL FLOW 


By 

G. Goglia^ and Sedat Biringen^ 

SUMMARY 

This paper involves a numerical simulation of the final stages of 
transition to turbulence in plane channel flow at a Reynolds number of 7500. 
Three^imensional, incompressible Navier-Stokes equations are numerically 
integrated to obtain the time-evolution of two- and three-dimensional fi- 
nite-amplitude disturbances. Computations are performed on the CYBER-203 
vector processor for a 32^33x32 grid. Solutions indicate the existence of 
structures similar to these observed in the laboratory and which are charac- 
teristic of various stages of transition that lead to final breakdown. 
Details of the resulting flow field after breakdown indicate the evolution 
of streak-like formations found in turbulent flows. Although the flow field 
does approach a steady-state (turbulent channel flow), implementation of 
subgrid-scale terms are necessary to obtain proper turbulent statistics. 

1. INTRODUCTION 

Recent experiments by Nishioka, Asia & lida (198l> have shown that 
transition to turbulence in a plane channel flow follows a sequence of 
events similar to that observed by Klebanoff, Tidstrom & Sargent (1962) in 
the boundary-layer transition. In this work, a direct numerical integration 
of the Navier-Stokes equations is performed in an attempt to simulate these 
events in plane channel flow, during the later stages of transition. 

In their experiments, Nishioka et al. (1981) used a vibrating ribbon 
technique to generate two-dimensional disturbances fixed at 72 Hz. To 

^Eminent Prof essor/Cha irman , Department of Mechanical Engineering and 
Mechanics, Old Dominion University, Norfolk, Virginia. 

. ... 

'“Research Associate Professor, Department of Mechanical Engineering and 

Mechanics, Old Dominion University, Norfolk, Virginia. 


excite the fully developed flow in a channel with a 27.4 aspect ratio. They 
measured the streamwise mean and fluctuating velocities, and u*, re- 

spectively, at a fixed streamwise location at a subcritical (linearly sta- 
ble) Reynolds number, Re * 5000 and simulated the various stages of transi- 
tion by varying the disturbance amplitude. Their observations show that 

subcritical instability takes place at a threshold amplitude of (u*) /U 

^ max 

* 0.01, where Uq is the mean velocity at the channel centerline. The 
evolution of this instability is evidenced by the intensification of the 
spanwise variation of the wavefront which develops into a peak-valley 
structure. Nishioka et al. (1981) observed that flow development follows 
then a trend which is similar to transition in the boundary-layer (Klebanoff 
et al. 1962, Kovasznay, Koraoda & Vasudeva 1962): local shear layers are 

t 

formed away from the wall at spanwise peak positions (ui/U^ • 0,11). In 


rapid succession, two-, three-, five- and multi-spike stages are observed 
with increasing amplitude of the primary disturbance. Nishioka et al. 

(1981) present evidence that in the final stages of transition, the flow 
starts to develop structures very similar to those found in fully developed 
wall turbulence. During this stage, the flow field is characterized by the 
development of a viscous sublayer, occurrence of the typical "streaks** close 
to the wall and the formation of horseshoe vortices sometimes referred to as 
the building blocks of wall turbulence (Theodorsen 1954). The present work 
attempts to simulate this sequence of events. 


Direct numerical integrations of the Navier-Stokes equations for the 
simulation of channel-flow transition has been the subject of some previous 
investigations. George & Heliums (1972) and Fasel, Bestek & Schefenacker 
(1977) used the two-dimensional Navier-Stokes equations to investigate the 
stability of channel flow to two-dimensional finite-amplitude disturbances. 
George & Helium (1972) studied the relationship between critical amplitude 
and Reynolds number and found a minimum Reynolds number, Re “ 3500, below 
which their predictions remained stable. This is contrary to experimental 
evidence (e.g., Kao & Park 1970) which indicates instability of plane chan- 
nel flow to finite-amplitude disturbances at Reynolds number as low as 1000 
Fasel et al. (1977) investigated the effects of disturbance amplitude on 
transition at subcritical and supercritical Reynolds numbers. They found 
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that increasing amplitude (uf/U^ = 0.06) of two-dimensional disturbances can 
drive plane channel flow to instability at a subcritical Reynolds number, Re 
» 5000. Although some insight into finite-amplitude instability of two- 
dimensional disturbances can be obtained from such calculations, a proper 
simulation of the transition process requires the u:)e of the three-dimen- 
sional Navier-Stokes equations. Once the process extends into the nonlinear 
regime, transition becomes increasingly three-dimensional so that for a 
relistic (physically plausible) representation of flow breakdown, spanwise 
variations of the flow field variables must be accounted for. It is also 
well known that creation of vorticity through vortex stretching, an essen- 
tial ingredient of transition, can take place only in a three-dimensional 
flow field. Hence, numerical solutions of the two-dimensional Navier-Stokes 
equations cannot represent energy transfer down the wave-number spectrum, 
which is the basic mechanisms of laminar flow transition to turbulence and 
the result of the vortex stretching mechanism. 

Effects of three dimensionality on transition have first been document- 
ed in detail by Klebanoff et al. (1962). Accordingly, three-dimensionality 
manifests itself mainly in the spanwise velocity variations resulting in the 
production of streamwise vorticity which, in turn, interacts with the span- 
wise vorticity and drives the flow to breakdown. Orszag & Kells (1980) and 
Patera & Orszag (1981) have expanded on this idea to study the susceptibil- 
ity of plane channel flow to three-dimensional disturbances by numerically 
integrating the three-dimensional Navier-Stokes equations. Hieir computa- 
tions at subcritical Reynolds numbers revealed some interesting aspects of 
subcritical transition. They found that initial disturbances, which are 
finite-amplitude two-dimensional Orr-Sommerfeld eigensolutions, decay slowly 
and, as expected, rate of decay increases with decreasing Reynolds nimiber. 
They also found that the addition of three-dimensional, finite-amplitude 
disturbances promote rapid instability at Reynolds numbers as low as 1000, 
which is in good agreement with the experiments of Kao & Park (1970). Their 
results suggest a similar tendency of the flow to instability even for 
small-amplitude, three-dimensional disturbances. They conclude that the 
mechanism that drives plane channel flow to instability is tie interaction 
of two-dimensional and three-dimensional disturbances, supporting the idea 
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that three-dimensionality is central to transition in plane channel flow. 

In a more recent work, Kleiser (1982) incorporated r spectral method to 
solve the three-dimensional Navier-Stokes equations starting with weakly 
three-dimensional initial conditions. Comparisons of his results with the 
experiments of Nishioka, Asai & lida (1980) up to the first spike stage are 
favorable even from a quantitative point of view supporting the idea that 
direct numerical simulations can be used as a mean of investigating the 
nonlinear transition process. 

It is apparent that the complete simulation of transition requires the 
inclusion of non-linearity and three-dimensionality as the fundamential 
characteristics of the flow field. These effects are of primary importance 
in the present work wherein the main emphasis is on the final stages of 
transition. Hence, for an accurate representation of the underlying physi- 
cal phenomena that take place during this stage of transition, the present 
simulation has been done by using the full three-dimensional, time-dependent 
Navier-Stokes equations. In order to drive the flow to instability and 
transition rapidly, calculations are performed at a linearly unstable 
Reynolds number (Re * 7500), with finite-amplitude two- and :hree-dimen- 
sional eigensolutions of the Orr-Soramerfeld equation used as the initial 
conditions. No attempt is made in this work to investigate the effects of 
different initial condition? or of Reynolds numbers. In section 2, the 
numerical methods used in the present study are briefly discussed. In 
section 3, results of calculations are presented and compared with the 
experiments of Nishioka et al. (1981). Finally, section 4 contains a sum- 
mary of results and some concludi*^g remarks. 


2. THE CALCULATION PROCEDURE 


The calculation procedure is based on the incompressible Navier-Stokes 
equations in primitive-variable form, 


9u . 
1 


+ 


9t 


9u . 
1 



- i i£_ + 



p ax. ax^ax^ 


( 1 ) 


and the continuity equation, 


4 


du . 

I 


0 


( 2 ) 


3x . 
1 


where are the velocities along the directions, p is the density, 

V is the kinematic viscosity and p is the total hydrostatic pressure. 

The equations are non‘-dimensionalized by the mean centerline velocity Uq 
and the channel half-width, h. The flow is assumed to be driven by a con- 
stant mean pressure gradient 2/Re, where Re is the Reynolds number given by 
Uoh/v. Also, the convective terms are put into a form which prevents 
occurrence of nonlinear instability in the numerical solution procedure by 
ensuring conservation of momentum and energy (Mansour, Ferziger, & Reynolds 
1978). The final form of the Navier-Stokes equations reads. 


3u- /9u. ^ 

1 fi i\ 3p2.^1 1 

+ j j ^ ^ — ^il 

9t ^^i/ ^'^i ^ Re 


(3) 


• I 

where P * p/p pressure head and 6. is the Kronecker 

delta. 


The flow is assumed to be periodic in the streamwise x^ and the span- 
wise x^ directions along which the flow field variables can be expanded in 
terms of Fourier series. This enables the use of the pseudo-spectral method 
(Orszag 1972) to calculate the spatial derivatives along x^ and x^ by 
use of discrete Fourier transforms. Considering transforms in the x 
direction, along which there are equally spaced mesh points, the veloc- 

ity component u^ can be written as 


Ui (x i) 


N/2-1 

^ ui(ki)e^ 

— N/2 


(4) 


where x^^ * mAxp m * 0,1,..., N-l and ■ 2IInjAX|. Accordingly, the Fourier 
transform of u^ is 


Ni-l 

u J ( k 1 ) * - — ^ Uj(xj)e 

N, m“0 


(5) 
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The spatial derivative of 


along can now be written as 


^1 


3ui(xi) 


dx 


1 


Ni/2 -1 

ikju^(kj)e^ 

ni*- Ni/2 


( 6 ) 


The derivative can be computed by forming the Fourier transform of u^(x^), 
multiplying the result by ik]^ and computing the inverse transform. For 
periodic functions, the pseudo-spectral method provides a means by which the 
spatial derivatives are evaluated with maximum accuracy for a given number 
of grid points. Along the X 2 *direction, a mesh stretching that concentrates 
grid points close to the solid walls is employed. The resulting mesh 
enables the resolution of the sublayer that is formed during transition for 
y+ < 2, where y^ is the coordinate along K 2 » in wall units. Spatial 
derivatives along X 2 are evaluated by a second-order finite-difference 
scheme on this non-uniform mesh. 

The governing equations were numerically integrated by the semi-implic- 
it method of Moin, Reynolds & Ferziger (1978), This procedure employs the 
explicit Adams-Bashforth method for the convective terms and the implicit 
Crank-Nicholson method for pressure and for the viscous diffusion terms. In 
order to start the two time-level Adams-Bashforth method, the Euler-implicit 
method is used at the first time step. 

Once the governing equations are discretized in time, a two-dimensional 
Fourier transform along the periodic directions x^ and X 3 transforms the 
equations into the kj-k 3 wave-number space. The transformed equations are 
written below in block-tridiagonal form for inversion along X 2 


A f"*; . B f"*‘ . C f"*; - r" 


(7) 


In (7), A, B, and £ are coefficient 
the advanced time level, n+ 1 , and at 
right-hand side vector that contains 
terms at the previous time levels. 


matrices, F. is the solution vector at 

n 

the X 2 “direc t ional node, j; Rj is the 
the convective, diffusive and pressure 
These are given as 
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Coefficients of the finite-difference operators that appear in the matrices 
A, and C are given as 

C2. - 2/A. Ja. , + A.). B2. - A2 . - 2/A.fA.^, + A.) 

J J + 1 J+J J J J*1 J J J J + 1 J 

Cl. - Al. - - A.) 

J J J 

“ ^* 2 ^ -...1 ~ (xo) • 

Since all the flow variables in the solution vector contain an imaginary and 
a real part, the block-inversion process is applied twice for each pair of 
and which are the wavenumbers along and x^, respectively. 

The assumption of periodicity in x^ and X 3 eliminates the necessi ty 
of applying explicit boundary conditions along these directions. However, 
due to the presence of solid boundaries along the X 2 “direction, no-slip 
boundary conditions are imposed on Up U 2 t ^nd U 3 and the pressure at 
the wall is calculated by a second order approximation from the interior of 
the flow field. That the pressure boundary conditions are consistent with 
the X 2 “momentum equation at the wall, has been shown by >toin et al. (1978). 

Initial conditions were prescribed from the two- and three-dimensional 
eigensolut ions of the Orr-Soramer feld equation by considering that even for 
subcritical Reynolds numbers, plane channel flow can be driven to instabili- 
ty if the least stable two-dimensional finite-amplitude Orr-Soramer fe Id 
eigenmodes are interacted with finite-amplitude three-dimensional eigenmodes 
(Orszag & Kells 1981). The most explosive situation arises when the three- 
dimensional eigenmodes are aligned with the main flow direction at ^45 to 
f60 degrees. Accordingly, we have used the following initial condition 

i Qtx I -i 0 x 3 

u(x) ■ U(x 2 , 0 , 0 ) + (X 2 )e u^^(x 2 )e^^^ ( 8 ) 
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Here, U(x2i0,0) is the parabolic velocity profile of plane channel flo/. 
The eigenfunctions U2 d(x 2) and ^ 3 d(x 2 ) correspond to two-dimensic .al and 

three-dimensional solutions of the Or r-Sommet feld equation at Re * 7500, 
respectively. The two-dimensional solution was obtained for a » 1 whereas, 
the three-dimensional solution was obtained for a«l;8-il. A computer 
program, which essentially uses the Kaplan filtering technique, was used for 
the solution of the Orr-Sommerf eld equation (Re3molds 1967). The final 
amplitudes were chosen so that the maximum value of the xj-directional two- 
dimensional disturbance was set equal to 0.1 IUq and the maximum amplitudes 
of the XI -directional three-dimensional disturbances were each set equal t^ 
O.OSUq. 


3. RESULTS AND DISCUSSION 


The finite-difference system (7) was solved on the CYBER-203 vector 
processor at NASA/Langley Research Center. A 32x33x32 mesh was employed 

along the xi x2 and X 3 -directions, respectively. The computer code was 
fully vectorized and vectorized library subroutines were used for the main 
computational operations that the solution teci'nique employs. These vector 
operations mainly are one-dimensional fast Fouri^^r transform (FFT) to cal- 
culate spatial derivatives with the pseudo-spectral method, two-dimensional 
FFT to transform the equations into "^3 wave-number space and block- 

tridiagonal matrix inversion along x2 . For the FFT operations, typical 
vector lengths were around 1000, which is an optimal vector length to take 
full advantage of the vector processor. For the block-tr idagonal matrix 
inversion (which is essentially a scalar operation) a vectorized subroutine 
that inverts a large number of tridiagonal systems simultaneously, was used. 
This procedure decreases CPU time significantly by reducing the number of 
scalar operations required to invert each system separately. The fully 
vectorized code takes about 5 sec. of CPU time per time step for the 32x33x 
32 mesh to solve the finite-difference system (7) on a computational box, in 
which the flow is confined between rigid walls at X£ ± 1. Periodicity 

lengths (box lengths) along xl and x3 were chosen so that the smallest 
wave numbers allowed in the computational domain were equal to a = 1 and 6 * 
1, respectively, i.e., the box length was set equal to I'n slong these 
d irect ions. 
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It should be recalled that the time-advancement scheme employed in this 
work is partly explicit (on the convective terms) and partly implicit (on 
the diffusion and pressure terms). Although in view of linear stabi lity 
analysis, implicit methods are unconditionally stable (extrapolation to 
nonlinear equations is sometime va*?ue), the mixed nature of the present 
scheme as well as the time-accurate nature of the problem under investiga- 
tion necessitate adherence to stability bounds of explicit schemes. There- 
fore, in all the calculations reported here, the convective stability condi- 
tion (the Courant-Fr iedr ichs-Lewy condition) that requires the Courant 
number (C.N.) to be always less than one and the diffusive stability condi- 
tion were obeyed. With (Ax^) . ■ 0.00694 and 0.02 < AT < 0.05, where AT 

is the non-dimensional time-step, through the course of the calculations to 
C.N. var ied as 


C.N. 



< 0.4 


(9) 


whereas the diffusive stability criterion, 


D, varied as 


D 


1 


Re 


AT 

(Ax)*- . 
min 


< 0.14 


( 10 ) 


so that the diffusive stability condition which requires D \ 0.5, was also 
always satisfied. 

In the subsequent parts of this section, results obtained from the 
numerical integration of the finite-difference system (7) for the evolution 
of the initial disturbances are compared with the experiments of Nishioka et 
al. (1981). It should, however, be noted that there are several differences 
existing between the conditions of this experiment and the present condi- 
tions. First, periodic boundary conditions employed in the computation 
along X| and x are not realized in the laboratory where the flow is 
periodic in time. Second, because of periodicity, the comput at iona 1 flow 
field evolves in time, not in space as in the laboratory. Oie jus t i f icat ion 
to the first difference can be advanced on the basis of previous numerical 
experiments, in which transition simulations of the flat-plate boundary 
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layer (Orszag 1976) and of the plane channel flow (Fasel et aK 1977) with 
proper inflow~out f low boundary conditions gave similar results to those 
simulations where periodic boundary conditions were applied. Therefore, it 
could be expected that the periodic boundary conditions of the present com- 
putations should not introduce any significant consequences for comparisons 
with the laboratory flow. The assumption of streamwise periodicity which 
implies the evolution of the flow in time (not in space) enables the most 
efficient use of available computer resolution by resolving only one wave- 
length. This assumption can be justified on the grounds that in an advected 
coordinate frame, the spatial evolution of the laboratory flo»’: is essential- 
ly equivalent to temporal growth (Orszag & Kells 1980). A third difference 
between the experiment and the computation is the Reynolds nimber. The 
laboratory flow has a subcritical Reynolds number. Re ■ 5000, whereas the 
computation was done at Re ■ 7500 which is linearly unstable for this flow. 
This selection of the Reynolds number was due to the necessity of forcing 
the computations into transition and breakdown with the least amount of 
computer expenses. It should be noted that in their experiments, Nishioka 
et al. obtained laminar channel flow up to Reynolds numbers around 9000, and 
found that wall phenomena which are characteristic of the final stages of 
transition are independent of Reynolds numbers. Hence, the difference 
between the Reynolds numbers of the experiment and the computation should 
not have any important consequences for the qualitative comparisons between 
the two sets of results. 

During the computations, data were stored at approximately every 100 
time-steps and results at these instances were compared with those of 
Nishioka et al. (1981) for the various stages of transition. The set of 
data obtained from the coiaputation that most closely matches a given stage 
in the laboratory flow is used for comparison with that stage in the experi- 
ment. Hence, comparisons are necessarily of a qualitative nature. 

In figure 1, a history of the time-evolution of the flow is given in 
terms of the maximum amplitude of the two-dimensional primary disturbance, 
its two-dimensional harmonic and the three-dimensional primary disturbance. 
The trends displayed by these quantities are generally similar to the 
I ults of Orszag & Kells (1980) which they obtained from computations per- 


formed at Re ■ 1250. The main features of these trends are the rapid 
decrease in the two-dimensional primary-wave amplitude and the rapid 
increase of the amplitude of its harmonic. Also» the three-dimensional 
primary-wave amplitude first increases, then decreases at around T ■ 40. In 
the present calculations, fluctuating variations of the amplitudes are 
observed to set in as early as T * 18 but, even at later times, the fluctua- 
tions are more controlled than those indicated by the results of Orszag & 
Kells (1980). That no ’’explosive** instabilities were found in the present 
computations is in accord with the findings of Nishioka, Asai & Tida (1980), 
which implies that breakdown in channel flow is as gradual as the growth of 
instabilities found in free shear flows. 

Plane-averaged (over the Xj^ - plane) flow-field quantities which are 
representative of various stages of transition are shown in figures 2-7. 
Plots of plane-averaged (mean) velocity profiles, Up are shown in figure 
2 for laminar (initial), late transition and ’’early-turbulence” stages at T 
» 0, T * 42 and T * 50, respectively. The u^ distribution at T * 50 has a 
strong resemblance to the turbulent channel flow profile, with increased 
velocity gradient at the wall and with a full profile indicative of turbu- 
lent mixing. Although, as expected, <u^> profiles do not show any fluctua- 
tions (or inflexions), plots of instantaneous velocity profiles do show very 
strong inflexions, especially in the regions close to the walls (figure 3). 
This indicates that the interaction of two- and three-dimensional waves 
close to the walls is the central mechanism that drives the flow to insta- 
bility. This is in accord with the idea that the flow will undergo transi- 

tion only for a selected band of spanwise wavenumbers, the most ’’dangerous” 
of which result in three-dimensional disturbances with maxima occurring 
close to the walls (Orszag & Kells I960). In figure 4, the velocity profile 

<u**“> * <u^>/u 2 versus y'*' X 2 U 2 /V is plotted; here U 2 is the friction 

velocity and is calculated from d <u^>/dx 2 at the wall. Although the plots 
indicate the formation of a sublayer, and the change from T * 42 to T * 50 
shows a gradual approach to the law-of-the-wall , the difference is still 
especially apparent in the logarithmic region. The values obtained from 
figure 8 of Nishioka et al. (1981) provides a fairly good comparison with 
the computation at T * 42 ; in fact, computational results at T * 42 will 
be presented as representing the five-spike stage of the experiment. At 
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At T ■ 42, the Reynolds number based on the friction velocity is equal to 
196 and, as expected, is considerably larger than its initial (laminar) 
value of 122 . 

Plots of plane-averaged fluctuation intensities, <(u^ - <u^>)2 >, are 
shown in figure 5 at various T* 'Hiere are several interesting features of 
this figure. Firstly, at T ■ 12, the shift in the position of peak ampli- 
tude towards the channel center (x 2 “ 0-4), as well as the increase in the 
maximum amplitude, indicate that the computation is developing in a manner 
compatible with the experimental observations pertaining to this stage of 
the transition process (Tani 1969). In addition to this, it will be shown 
later that there is a substantial increase in the spanwise vorticity, 
away from the wall at around X 2 ® 0.4. Secondly, in accordance with the 
laboratory flow of Nishioka et al. (1981) at later stages in computation 
(e.g. at T * 42), the intensity profiles display a second peak occurring 
close to the wall which sho\ild be associated with turbulence production. 
However, in spite of these similarities, a comparison of the computational 
results at T * 42 and the five-spike stage of the experiment (Nishioka et 
al.) shows that the locations of the peak positions do not coincide and the 
peak value obtained in the experiment is about 15 percent greater than that 
obtained from the computation. 

Spanwise variations of U]^ and 03 at T “ 42 are plotted at various 
distances along in figures 6 and 7, respectively. These figures 

indicate that at this stage the flow is highly three-dimensional; however, 
the spanwise symmetry imposed by the initial conditions is slill retained. 

An estimate of the flow field resolution along X 3 can be obtained from the 

spanwise distance between peak position, X. Nond imens ional ized by u2 find 
V typically X = 120. This is larger than but comparable to X =» 80, which 
is the typical spanwise length in the laboratory flow during the five— spike 
stage (Nishioka et al. 1981). It should be noted that the spanwise charac- 
teristic length in wall turbulence is about 100. It could, therefore, be 
asserted that, at this stage present results are representative of initial 
wall turbulence, 

A more detailed description of the transition process can be obtained 

9u I 

from contour plots of equi-shear lines, (which correspond to approximate 

3x2 
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spanwise vorticity, in Che - X 2 plane at the position of maximum 

ui/Uq. Computational result? that correspond to the various stages of the 
laboratory flow are presented in figures 8-12, In figures 9a * 11a, figures 
4-6 of Nishioka et al. (1981) are also presented for comparison with the 
present computations. In figure 8, contour plots corresponding to the ini- 
tial conditions and in figure 9, contour plots corresponding to the "one 
spike" stage (T ■ 0 and T ■ 12, nispectively) are sho%m. In ooth the 
laboratory flow and the computation, the typical head of the shear layer 
appears very clearly at the one-spike stage, indicating the formation of a 
shear layer away from Che wall at about X£ • 0.4 due to the induced velocity 
from the strearawise vortex system. In the experiment, the sudden dip of the 
shear layer from Che high-velocity outer flow to the low velocity region 
appears as a kink. This is not very clear in the computation and can be 
attributed to inadequate resolution along X 2 in the outer portions (close 
to the centerline) of the flow field. Inadequate mesh spacing in these 
portions of the flow field away from the wall manifests itself in the com- 
parison of contour levels and concentrations of approximate spanwise 
vorticity inside the "head" of the shear layer. This comparison mainly 
reveals that the laboratory flow exhibits higher vorticity levels than are 
indicated by the computation. However, since the mesh is finely clustered 
aloiig X 2 close to the wall, vorticity concentrations in this region are 
adequately resolved by the numerical simulation. 

Figure 10 shows equi-shear lines at T * 22 corresponding to the three- 
spike stage of the laboratory flow. In both the experiment and the computa- 
tion, due to the secondary instability manifested in the previous stage, 
breakdown of fl3w structures into smaller scales are observed. In the 
experiment, the growth of the kinked portions of the equi-shear lines into 
the so-called "hairpin eddies" is clearly depicted. The computation dis- 
plays a similar evolution, with the head of the shear layer lifted up to- 
wards the channel centerline, and the kink in the shear layer being quite 
apparent in this stage of the computation. Simultaneously with this activi- 
ty taking place in the outer (high speed) portions of the flow field, both 
the experiment and the computation show an intense shear layer developing 
close to the wall, which is indicative of turbulence generation. Although 
it is generally agreed that hairpin eddies which are lifted towards the 
centerline erupt into turbulent spots, it seems likely that wall-turbulence 
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is also closely associated with vorticity dynamics that take place simulta- 
neously with the eruption of hairpin eddies. Contours of equi-shear lines 
at T ■ 42 corresponding to the five-spike stage of the experiment are shown 
in figure 11. The intense shear-layer developed in the wall region is ap- 
parent in both the experiment and the computation. Vorticity levels in the 
laboratory flow and the computation are comparable but the laboratory flow 
indicates generally higher levels. The most significant feature in figure 
11 is the existence of distinct vortex structures in the wall region both in 
the laboratory flow and in the computation. These vortices are aligned 
close to 45* to the mean flow direction and show a close resemblance to the 
energetic horseshoe vortices which are characteristic of wall-turbulence. 

It is these vortices (turbulent eddies) that are mainly responsible for 
extracting energy from the mean shear (Tennekes & Lumley 1974, p, 41). In 
figure 12, equi-shear lines are shown at T “ 50. Here, all high-shear con- 
centration and intense vorticity dynamics are confined to the wall region. 
Also, :»w field is characterized by scales which are much smaller than 

those observed in the earlier stages, somewhat indicating that the flow 
field at this stage assumes a turbulent-like character. 

Contours of streamwise vorticity, o)x, are shown in figures 13-17 at 
various instances in time in the x2 - x3 plane at the position of maximum 
^ Figure 13 displays that at T =* 0, the counter-rotating vortex 

system which is due to the prescribed initial conditions is relatively weak 
and there are no vorticity concentrations close to the wall. Figure 14 

shows ojx at T » 12, corresponding to the one-spike stage of the labora- 
tory flow. At this stage in the simulation, the initial counter-rotating 
vortex system is still evident and there is an enhancement and concentration 
of streamwise vorticity close to the wall. Streamwise vorticity contours 
corresponding to the three-spike stage of the laboratory flow are shown in 
figures 15 at T * 22, where £g enhanced and its magnitude is compara- 
ble to the lateral component, In addition to this enhancement of 

0)x, figure 15 displays occurrence of patches of vorticity of scales much 
smaller than the initial counter-rotating vortices. Formation of alternat- 
ing positive (solid lines) and negative (dashed lines) vorticity regions 
close to the wall is also evident from this figure. Figure 16 and 17 show 
contours of streamwise vorticity during the multi-spike and early-turbulence 
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stages^ i.e. at T * 42 and T * 50, respectively. In these figures, the 
major features are further enhancement of streainwise vorticity, breakdown 
into smaller scales and concentration of streainwise vorticity close to the 
walls. Details of contours in the vicinity of the wall are shown in 

figure 18 at T “ 50 where the region between the lower wall (x2 ■ -1.0) and 
X 2 « -0.92 is magnified. The most interesting aspect of this figure is the 
distinct pattern of alternating regions of very high concentrations of posi- 
tive and very high concentrations of negative vorticity; however, due to 
magnification contour lines are distorted in this figure. The alternating 
streamwise-vortici ty region is usually associated with an induced velocity 
field that is responsible for vorticity production during the final stages 
of the transition process (Komoda 1967). It should also be noted that the 
alternating structure of positive and negative concentrations of vorticity 
close to the wall is a very distinct feature of fully developed wall-turbu- 
lence (Moin & Kim 1981). 

The mechanism that is responsible for the generation of vorticity con- 
centrations is usually explained as vortex stretching (and deformation) by 
the mean flow. The stretching and deformation moves downstream with a 
translation velocity that induces lower local mean velocity of the upstream 
edge of the vorticity layer than its downstream edge (Komoda 1967). An 
examination of contours of mean velocity (figure 19) along with the spanwise 
vorticity contours at T » 22 (figure 10) in the Xj - X 2 plane at the posi- 
tion of maximum /Uq clearly shows a similar trend. The nose of the 
negative streamwise vorticity layer (which corresponds approximately to 
regions with high concentrations of equi-shear lines) is generally associ- 
ated with higher velocities than the upstream region and large variations of 
the local mean velocity exist within the layer. 

Normal velocity, U 2 , contours are shown in figures 21-24 during the 
three-spike, multi-spike and ear ly-turbulence stages at the position of 
maximum u^ /Uq. Figure^ 20 shows U 2 -contours in the x\ - X 2 plane at T * 

22, corresponding to the three-spike stage of the laboratory flow. In this 
figure, there is clear evidence of the beginning of an alternating up- and 
-down flow similar to the pattern described by Kovasznay et al. (1967), 
which is characterized by the intense updrift accompanied by fluid drifting 
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down at both sides. This alternating structure of up-and-down flow in the 
Xj - X 2 plane is more evident at T « 42 (figure 21), corresponding to the 
multi-spike stage of the experiment. Also in this figure, alternating 
regions of fluid with scales much smaller than those at T * 22, are clearly 
depicted. In figure 22 and 23, normal velocity contours are shown in the 
X 2 • X 3 plane at T * 42 and T ■ 50 corresponding to the multi-spike and 
early-turbulence stages, respectively. These contours are plotted in the 
region between the lower wall (x 2 * -1.0) and X 2 • *0.92. In these figures, 
the development of the up-and-down pattern of the fluid close to the wall is 
observed; the contours at T * 50 display alternating structures very similar 
to the characteristic streak-like structures found in the wall region of 
turbulent channel flow (Mo in & Kim 1981). 

4. SUMMARY AND CONCLUDING REMARKS 

In this study, final stages of transition to turbulence in plane chan- 
nel flow have been simulated by a direct numerical solution of the Navier- 
Stokes equations. Results were compared with the experiments of Nishioka et 
al. (1981) and it was found that, in spite of the limited resolution of the 
32x33x32 grid employed in the computations, the simulation was capable of 
reproducing most of the essential features of wall phenomena observed in the 
laboratory flow. Grid resolution in the Xj- and x 3 -directions, along v’nich 
the flow is periodic, was found to be adequate to capture the sequence of 
events that lead to early turbulence. Vorticity and velocity contours in 
the vicinity of the lower wall indicated formation of horseshoe vortices 
and, at later stages, streak-like structures alternating in the spanwise 
direction. Typically, the spanwise characteristic length, X, inferred 
from the spanwise variations of Uj and U 3 , was found to be around X ■ 
120, which is close to X * 100 of fully developed wall turbu'ance. In 
agreement with the experiment, it was found that during the later stages of 
transition, flow field statistics indicate the formation of a laminar 
sublayer; however, the development of the logarithmic region and hence the 
approach to fully developed turbulence is slow. 

The main deficiency of this study inevitably stems from limited spatial 
resolution and manifests itself in several ways. Firstly, at later stages 
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of the computation, e.g. for T ■ 42, the computation was not able to main- 
tain the total disturbance energy (turbulent kinetic energy). This is not 
surprising and can be explained on the grounds that close to the wall, the 
available mesh is not able to resolve the characteristic large-scale struc- 
tures, which appear as finely spaced streaks. Insufficient mesh resolution 
could also result in lower gradients of the mean velocity in the viscous 
sublayer. The combined effect of these is less turbulence production and 
therefore indefinitely decaying turbulent kinetic energy. Secondly, the 
finite cut-off wave numbers along xx and X 3 prevent energy transfer down 
the wave-number space. Accordingly, at large T, pollution of Fourier 
modes and accumulation of excess energy at low wave numbers become very 
significant sources of inaccuracy. Hence, proper simulation of transition 
beyond the early turbulence stage necessitates the use of the higher grid 
resolution and even then, the incorporation of subgrid scale turbulence 

modeling for an accurate and realistic representation of the flow field 
phenomena. 
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Reynolds for helpful discussion and to P.J. Bobbitt and W.D. Harvey for 
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Figure 6. Spanwiae variations of Uj at T = 42. (a) X 2 = 0.029; (b) 

0.099; (c) X, = 0.256. 
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Figure 10. Contour plots of - — . (a) three-spike stage, figure from 

3x2 

Nishioka et al. (1981); (b) computations at T = 22 in the Xj,X 2 
plane; contours from -0.3 to 5.1. 


30 













ORIGINAL PAGE IS 
OF POOR QUALITY 



9ui 

Figure 12. Contour plots of , at T * 50 in the x^X 2 plane; contours 

9x2 

from -0.6 to 9.0. 
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Figure 14 * Contour plots of in the x-. X3 plane at T * 12 ; contours 

from -1*6 to 1.6. 
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Figure 15. Contour plots of 
from “4.8 to 4 . 8 . 


in the ^ 2 * ^3 plsne at T ■ 22; contours 
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Figure 19. 
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Contour plots of 
from -0. to 0.96. 


in the Xp X 2 plane at T 


22 ; contours 
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Figure 20. Contour plots of U 2 in the Xj, Xj plane at T » 22; contours 
from -0.081 to 0.072 labels scaled by 1000. 
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Figure 23. 


Contour plots of U£ in the X 2 , Xj 
vicinity of the lower wall; contours 


plane at T ■ 50 
from -0.135 to 0 


in 
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